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Abstract 

The graphical lasso [Friedman et al.l . 120071 ] is an algorithm for learning the structure in an 
undirected Gaussian graphic al model, using i\ regularization to control t he number of zeros in the 
precision matrix = S _1 [Baneriee et all 120081 . lYuan and Linl . l2007f | . The R package GLASSO 
[Friedman et al.l . |2007( | is popular, fast, and allows one to efficiently build a path of models for 
different values of the tuning parameter. Convergence of GLASSO can be tricky; the converged 
precision matrix might not be the inverse of the estimated covariance, and occasionally it fails to 
converge with warm starts. In this paper we explain this behavior, and propose new algorithms 
that appear to outperform GLASSO. 

By studying the "normal equations" we see that, GLASSO is solving the dual of the graphi- 
cal lasso penal i zed li kelihood, by block coordinate ascent; a result which can also be found in 
iBaneriee et al.l [20081 ] . In this dual, the target of estimation is S, the covariance matrix, rather 
than the precision matrix 0. We propose similar primal algorithms P-GLASSO and dp-GLASSO, 
that also operate by block-coordinate descent, where is the optimization target. We study all of 
these algorithms, and in particular different approaches to solving their coordinate sub-problems. 
We conclude that dp-GLASSO is superior from several points of view. 



1 Introduction 



Consider a data matrix X nxp , a sample of n realizations from a p-dimensional Gaussian distribution 
with zero mean and positive definite covariance matrix S. The task is to estimate the unknown S 
based on the n samples — a challenging problem especially when n <C p, when the ordinary maximum 
likelihood estimate does not exist. Even if it does exist (fo r p < n), the MLE is of ten poorly behaved, 
and regularization is called for. The Graphical Lasso [Friedman et all 12007 is a regularization 
framewor k for estimating the covariance matrix S, under the assumption that it s inver se = E _1 



is sparse Baneriee et al. , 20081 Yuan and Lin . 2007 . Meinshausen and Biihlmann . 2006| . is called 
the precision matrix; if an element 9jk = 0, this implies that the corresponding variables Xj and Xk 
are conditionally independent, given the rest. Our algorithms focus either on the restricted version 
of or its inverse W = _1 . The graphical lasso problem minimizes a ^-regularized negative 
log-likelihood: 

minimize/(0) := -logdct(0) + tr(S0) + A[|0[|i. (1) 

Here S is the sample covariance matrix, ||0||i denotes the sum of the absolute values of 0, and A 
is a tuning parameter controlling the amount of £-\ shrinkage. T his is a semidefinite programming 



problem (SDP) in the variable [Bovd and Vandenberghd . 12004 



In this paper we revisit the GLASSO algorithm proposed by iFriedman et al.l [2007| for solving (p}; 
we analyze its properties, expose problems and issues, and propose alternative algorithms more 
suitable for the task. 
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1 



Some of the results and conclusions of this paper can be found in iBaneriee et al.l [20081 ] , both 
explicitly and implicitly. We re-derive some of the results and derive new results, insights and 
algorithms, using a unified and more elementary framework. 



Notation We denote the entries of a matrix A nxn by ay. ||A||i denotes the sum of its absolute 
values, ||-A||oo the maximum absolute value of its entries, ||^4||_f is its Frobcnius norm, and abs(A) is 
the matrix with elements |a^ |. For a vector u G 3? 9 , ||u||i denotes the l\ norm, and so on. 
From now on, unless otherwise specified, we will assume that A > 0. 



2 Review of the GLASSO algorithm. 



We use the frame- work of "normal equations" as in lHastie et ail |2009j | . iFriedman et al.l [2007} . Using 



sub-gradient notation, we can write the optimality conditions (aka "normal equations" ) for a solution 
to (JTJ as 

-0 _1 + S + Ar = O, (2) 
where T is a matrix of component- wise signs of 0: 

"f jk = sign(^ fe ) if 9 jk ^ , > 

ljk G [— 1, 1] if ©j-fc = ( > 

(we use the notation 7^ G Sign(flj-fc)). Since the global stationary conditions of ([2]) require 9jj to be 
positive, this implies that 

wu = su + A, % = 1, . . . ,p, (4) 

where W = _ . 

GLASSO uses a block-coordinate method for solving ([2]). Consider a partitioning of and T: 

0ii 012 \ r _f r u 7ia 



e = « « . r = ( 5 ) 

V 021 ^22 / V 721 722 / 

where 0n is (p — 1) x (p — 1), #12 is (p — 1) x 1 and #22 is scalar. W and S are partitioned the 
same way. Using properties of inverses of block-partitioned matrices, observe that W = 1 can be 
written in two equivalent forms: 



Wn W12 

W21 W22 I 



012021 W W 012 

(6) 



Q-l _|_ ©11 012^21©n ®n #1 



#22— 021©n ^12) ^22— ^21©n #12 



?22— ^21©n 1 ^12) 



(7) 



GLASSO solves for a row/column of ([2]) at a time, holding the rest fixed. Considering the pth column 
of (J2J, we get 

- W12 + S12 + A7 12 = 0. (8) 

Reading off W12 from ((6]) we have 

W12 = -Wll0l2/022 (9) 

and plugging into ©, we have: 

W 11 ^+s 12 + A 7l2 = 0. (10) 

P22 

GLASSO operates on the above gradient equation, as described below. 
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As a variation consider reading off W12 from iJTj): 



Sl2 + A 7 12=0. (11) 



(#22 — ^21©n 1 ^12) 

The above simplifies to 

n 1 0i2W22 + si2 + A7 12 =0, (12) 

where w 22 = 1/(6*22 — 02i©i~i ^12) i s fixed (by the global stationary conditions (TJ|). We will see that 
these two apparently similar estimating equations (|10[) and (|12[) lead to very different algorithms. 
The GLASSO algorithm solves ([TO]) for (3 = 9 12 /9 22 , that is 

Wii/3 + Si 2 + A 7l2 =0, (13) 

where -f 12 £ Sign(/3), since 6* 2 2 > 0. (fT"3"|) is the stationarity equation for the following £1 regularized 
quadratic program: 

minimize {±/3'Wn/3 + /3'si 2 + A||/3||i} , (14) 

/3S5RP- 1 

where Wn >- is assumed to be fixed. This is analogous to a lasso regression problem of the last 
variable on the rest, except the cross-product matrix Sn is replaced by its current estimate Wn. 
This problem itself can be solved efficiently using elementwise coordinate descent, exploiting the 
sparsity in (3. From /3, it is easy to obtain W12 from ©. Using the lower-right element of ©, 6*22 is 
obtained by 

= W22 — ft W12. (15) 

#22 

Finally, 12 can now be recovered from and 6*22 ■ Notice, however, that having solved for (3 and 
updated W12, GLASSO can move onto the next block; disentangling 9\ 2 and #22 can be done at the end, 
when the algorithm over all blocks has converged. The GLASSO algorithm is outlined in Algorithm Q] 
We show in Lemma [3] in Section |8] that the successive updates in GLASSO keep W positive definite. 



Algorithm 1 GLASSO algorithm [Friedman et all 12007 



1. Initialize W = S + AI. 

2. Cycle around the columns repeatedly, performing the following steps till convergence: 

(a) Rearrange the rows/columns so that the target column is last (implicitly). 

(b) Solve the lasso problem (|14p . using as warm starts the solution from the previous round 
for this column. 

(c) Update the row/column (off-diagonal) of the covariance using W12 ©. 

(d) Save f3 for this column in the matrix B. 

3. Finally, for every row/column, compute the diagonal entries 9jj using (|15p . and convert the B 
matrix to 0. 



Figure[T](left panel, black curve) plots the objective /(© ) for the sequence of solutions produced 
by GLASSO on an example. Surprisingly, the curve is not monotone decreasing, as confirmed by the 
middle plot. If GLASSO were solving (JlJ by block coordinate-descent, we would not anticipate this 
behavior. 

A closer look at steps (j9j) and (fT0|) of the GLASSO algorithm leads to the following observations: 

(a) We wish to solve ([8]) for 9i 2 . However 9i 2 is entangled in Wu, which is (incorrectly) treated 
as a constant. 

(b) After updating 9i 2 , we see from ([7]) that the entire (working) covariance matrix W changes. 
GLASSO however updates only W12 and W21. 
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Figure 1: [Left panel] The objective values of the primal criterion (Jj) and the dual criterion I19(l corre- 
sponding to the covariance matrix W produced by GLASSO algorithm as a function of the iteration index ( each 
column/row update). [Middle Panel] The successive differences of the primal objective values — the zero 
crossings indicate non-monotonicity. [Right Panel] The successive differences in the dual objective values — 
there are no zero crossings, indicating that GLASSO produces a monotone sequence of dual objective values. 



These two observations explain the non-monotone behavior of GLASSO in minimizing /(©). Scction[3] 
shows a corrected block-coordinate descent algorithm for 0, and Section 2] shows that the GLASSO 
algorithm is actually optimizing the dual of problem (TTJ), with the optimization variable being W. 



3 A Corrected GLASSO block coordinate-descent algorithm 

Recall that ([T2jl is a variant of (fT0|) . where the dependence of the covariance sub-matrix Wn on 9\2 
is explicit. With a. = 6J12W22 (with W22 > fixed), 811 >- 0, (fl~2|) is equivalent to the stationary 
condition for 

minimize {ia'07/a + a'si9 + A||a||i) . (16) 

aesfp- 1 12 

If a is the minimizcr of (|16|) . then 612 = a/u'22- To complete the optimization for the entire 
row/column we need to update 622- This follows simply from ([7]) 

022 = — + 021©n 012, (17) 

W22 

with W22 = S22 + A. 

To solve (fT6|) we need &J{ for each block update. We achieve this by maintaining W = _1 as 
the iterations proceed. Then for each block 

• we obtain &11 from 

©u = w n _ W12W21/W22; (18) 

• once 612 is updated, the entire working covariance matrix W is updated (in particular the 
portions Wn and W12), via the identities in ([7]), using the known 

Both these steps are simple rank-one updates with a total cost of 0(p 2 ) operations. 

We refer to this as the primal graphical lasso or P-GLASSO, which we present in Algorithm [2 
The P-GLASSO algorithm requires slightly more work than GLASSO, since an additional 0(p 2 ) 
operations have to be performed before and after each block update. In return we have that after 
every row/column update, and W are positive definite (for A > 0) and 0W = I p . 
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Algorithm 2 P-GLASSO Algorithm 



1. Initialize W = diag(S) + AI, and = W" 1 . 

2. Cycle around the columns repeatedly, performing the following steps till convergence: 

(a) Rearrange the rows/columns so that the target column is last (implicitly). 

(b) Compute 0^ using (T8|). 

(c) Solve (|16p for a, using as warm starts the solution from the previous round of row/column 
updates. Update 6i 2 = &/W22, and 622 using (fTT)) . 

(d) Update and W using ([7]), ensuring that 0W = I p . 

3. Output the solution (precision) and its exact inverse W (covariance). 



4 What is GLASSO actually solving? 

Building upon the framework developed in Section [21 we now proceed to establish that GLASSO 
solves the convex dual of problem ((T|), by block coordinate ascent. We reach this conclusion via 
elementary arguments, closely aligned with the framework we develop in Section [5J The approach 
we present here is intended for an audience without much of a familiarity with convex duality the- 
ory |B^yd_aird_Vanderiber^heJ 2004 1. 



Figure [T] illustrates that GLASSO is an ascent algorithm on the dual of the problem [TJ The red 
curve in the left plot shows the dual objective rising monotoncly, and the rightmost plot shows that 
the increments are indeed positive. There is an added twist though: in solving the block-coordinate 
update, GLASSO solves instead the dual of that subproblem. 

4.1 Dual of the l\ regularized log-likelihood 



We present below the following lemma, the conclusion of which also appears in lBaneriee et all [20081 ] . 
but we use the framework developed in Section [2j 

Lemma 1. Consider the primal problem fl]) and its stationarity conditions These are equivalent 
to the stationarity conditions for the box- constrained SDP 

maximize 5(f) := log det(S + f ) + p (19) 

f: ||f ||oo<A 

under the transformation S + T = ®~ . 

Proof. The (sub)gradient conditions ((2j) can be rewritten as: 

- (S + AT)" 1 +0 = (20) 

where T = sgn(0). We write T = XT and observe that ||r||oo < A. Denote by abs(0) the matrix 
with element- wise absolute values. 

Hence if (0,T) satisfy (f20"|) . the substitutions 



T = AT; P = abs(0) (21) 

satisfy the following set of equations: 

-(S + fT 1 +P*sgn(f) = 

P * (abs(f ) - AlglJ,) = (22) 
l|f||oo < A. 

In the above, P is a symmetric pxp matrix with non- negative entries, lplp denotes a pxp matrix of 
ones, and the operator '*' denotes element-wise product. We observe that ([2^]) arc the KKT optimality 
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conditions for the box-constrained SDP (|T9")) . Similarly, the transformations = P * sgn(f) and 
r = r/A show that conditions (|2"2"j) imply condition (|2T))) . Based on (|20[) the optimal solutions of the 
two problems (P and JTSJ) are related by S + f = _1 . □ 

Notice that for the dual, the optimization variable is T, with S + T = _1 = W. In other words, 
the dual problem solves for W rather than 0, a fact that is suggested by the GLASSO algorithm. 

Remark 1. The equivalence of the solutions to problems U9\) an d ([7]) as described above can also 
be derived via convex duality theory \Boud and Vandenberahe], 200A l. which shows that U9\) is a dual 



function of the l\ regularized negativ e log-likelihood 



solutions of the two problems coincide \Baneriee et al 



Strong duality holds, hence the optimal 
'200& 1. 



We now consider solving for the last block j 12 (excluding diagonal), holding the rest of f 
fixed. The corresponding equations are 

-012 + P12 * sgn(7 12 ) = 
P12 * (abs(7 12 ) - Alp-i) = (23) 
Il7i 2 lloo < A. 

The only non-trivial translation is the 612 in the first equation. We must express this in terms of the 
optimization variable 7 12 - Since Si2+7 12 = W12, using the identities in ©, we have W7/] 1 (si2+7 12 ) = 
— 612/622- Since 822 > 0, we can redefine pi 2 = Pi2/#22, to get 

Wi 1 1 (si 2 + 7 12 ) + P12 * sgn(7 12 ) = 

P12 * (abs(7 12 ) - Mp-i) = (24) 
II712II00 < A. 

The following lemma shows that a block update of GLASSO solves (jM]) (and hence (J23J)), a block 
of stationary conditions for the dual of the graphical lasso problem. Curiously, GLASSO does this not 
directly, but by solving the dual of the QP corresponding to this block of equations. 

Lemma 2. Assume Wn >- 0. The stationarity equations 

Wiijl + si2 + A 7l2 = 0, (25) 
where 7 12 € Sign(f3), correspond to the solution of the l\-regularized QP: 

minimize ±/3'W n /3 + /3'si 2 + A||/3||i. (26) 

Solving V2b}) is eguivalent to solving the following box- constrained QP: 

minimize i(s 12 + 7)'W n 1 (s 12 + 7) subject to ||'y||oo < A, (27) 



with stationarity conditions given by {24]) , where the (3 and 7 12 are related by 

^ = -W n 1 (s 12 +7i2)- (28) 
Proof. (|25|) is the KKT optimality condition for the l\ regularized QP (|26| . We rewrite (|25|) as 

^ + W n 1 ( Sl 2 + A7 12 ) = 0. (29) 

Observe that $i — sgn(/3j)|/3j| Vi and ||7 12 ||oo < 1. Suppose f3, 7 12 satisfy (|2"9"|) . then the substitutions 

7i2 = A7i2> P12 = abs0) (30) 

in ([2"9"| satisfy the stationarity conditions (f2"4"]). It turns out that §M§ is equivalent to the KKT 
optimality conditions of the box-constrained QP (|27| . Similarly, we note that if 7i 2 ,Pi2 satisfy ([24]) , 
then the substitution 

7i2 = 7i 2 / A 5 h = P12 * sgn(7 12 ) 
satisfies ([2^1) . Hence the and 7 12 are related by (|28l) . □ 
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Rem ark 2. The above result can also be derived via convex duality theorv tBoud and Vandenberaht , 
200A l. where \21ty is actually the Lagrange dual of the l\ regularized QP V2b}). with \28\) denot- 
ing the primal-dual relationship. \Baneriee et al\ . \2004 . Section 3.3] interpret \2T§ as an l\ pe- 
n alized regression problem (using co nvex duality theory) and explore connections with the set up 
of \Meinshausen and BuhlmanA \200d l. 

Note that the QP ([27]) is a (partial) optimization over the variable W12 only (since S12 is fixed); 
the sub- matrix Wn remains fixed in the QP. Exactly one row/column of W changes when the block- 
coordinate algorithm of GLASSO moves to a new row/column, unlike an explicit full matrix update 
in Wn, which is required if #12 is updated. This again emphasizes that GLASSO is operating on the 
covariance matrix instead of 0. We thus arrive at the following conclusion: 



Theorem 1. GLASSO performs block- coordinate ascent on the box- constrained SDP M9j) . the Lagrange 
dual of the primal problem |ip. Each of the block steps are themselves box- constrained QPs, which 
GLASSO optimizes via their Lagrange duals. 

In our annotation perhaps GLASSO should be ca lled DD-GLASSO. sinc e it performs dual block 
updates for the dual of the graphical lasso problem . iBaneriee et al.l [2008| , the paper that inspired 
the original GLASSO article [Friedman et all 12007 1. also operates on the dual. They however solve 



the block-updates directly (which are box constrained QPs) using interior-point methods. 

5 A New Algorithm — DP-GLASSO 

In Section [31 we described P-GLASSO, a primal coordinate-descent method. For every row/column we 
need to solve a lasso problem (|16[) . which operates on a quadratic form corresponding to the square 
matrix &7{". There are two problems with this approach: 

• the matrix &7-1 needs to be constructed at every row/column update with complexity 0{p 2 ); 

• &7i is dense. 

We now show how a simple modification of the £i-rcgularizcd QP leads to a box-constrained QP with 
attractive computational properties. 

The KKT optimality conditions for ([T5|) . following (fT2")l . can be written as: 

e^a + S12 + Asgn(a) = 0. (31) 

Along the same lines of the derivations used in Lemma [U the condition above is equivalent to 

qi2 *sgn(7) + ©ii (S12 + 7) = 

q X2 * (abs( 7 ) - Alp_i) = (32) 

Il7lloo < A 

for some vector (with non- negative entries) qi2- (|32[) arc the KKT optimality conditions for the 
following box-constrained QP: 

minimize ±(si 2 + 7)'0ii(si 2 + 7); subject to H7H00 < A. (33) 



The optimal solutions of ([33)) and (j31[) are related by 

a = -©11(812+7), (34) 

a consequence of (j31[) . with a = 612 ■ W22 and W22 = S22 + A. The diagonal #22 of the precision matrix 
is updated via (J7J): 

22 = 1 - (S12+ ^ 12 (35) 
W22 
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Algorithm 3 dp-GLASSO algorithm 



1. Initialize = diag(S + AI) _1 . 

2. Cycle around the columns repeatedly, performing the following steps till convergence: 

(a) Rearrange the rows/columns so that the target column is last (implicitly). 

(b) Solve ([53]) for 7 and update 

O12 = -011 (S12 + j)/w 2 2 

(c) Solve for 622 using (|35|) . 

(d) Update the working covariance W12 = S12 + 7. 



By strong duality, the box-constrained QP (|3"3"[) with its optimality conditions ([3"2"[) is equivalent 
to the lasso problem (|16[) . Now both the problems listed at the beginning of the section are removed. 
The problem matrix ©n is sparse, and no 0(p 2 ) updating is required after each block. 

The solutions returned at step 2(b) for 612 need not be exactly sparse, even though it purports 
to produce the solution to the primal block problem (|16[) . which is sparse. One needs to use a tight 
convergence criterion when solving (|33|) . In addition, one can threshold those elements of 812 for 
which 7 is away from the box boundary, since those values are known to be zero. 

Note that dp-glasso does to the primal formulation (J]) what glasso does to the dual, dp- 
GLASSO operates on the precision matrix, whereas GLASSO operates on the covariance matrix. 



6 Computational Costs in Solving the Block QPs 

The i\ regularized QPs appearing in (fT4|) and (fT6| are of the generic form 



minimize |u'Au + a'u + Allulli, (36) 

for A >- 0. In this paper, we choose to use c yclical coordinat e desc ent for solving (f36|) . as it is used 
in the GLASSO algorithm implementation of iFriedman et aL [2007 1 . Moreover, cyclical coordinate 



descent methods perform well with good warm-starts. These are available for both fTJJ and (pi 
since they both maintain working copies of the precision matrix, updated after every row/column 
update. There are other efficient ways f or solving (|3"6"1). capable of scaling to lar ge problems - 
for example first-order proximal methods (Beck and Tebouilel 20091 Nesterov . 2007 . but we do not 
pursue them in this paper. 

The box-constrained QPs appearing in ([27]) and (|3"3"[) are of the generic form: 

minimize |(v + b)'A(v + b) subject to HvH^ < A (37) 

for some A y 0. As in the case above, we will use cyclical coordinate-descent for optimizing ([371) . 

In general it is more efficient to solve ([3"6"[) than ([37) for larger values of A. This is because a large 
value of A in (pJo) results in sparse solutions u; the coordinate descent algorithm can easily detect 
when a zero stays zero, and no further work gets done for that coordinate on that pass. If the solution 
to ([36[) has « non-zeros, then on average « coordinates need to be updated. This leads to a cost of 
0(qn), for one full sweep across all the q coordinates. 

On the other hand, a large A for ([37]) corresponds to a weakly-regularized solution. Cyclical 
coordinate procedures for this task are not as effective. Every coordinate update of v results in 
updating the gradient, which requires adding a scalar multiple of a column of A. If A is dense, this 
leads to a cost of 0(q), and for one full cycle across all the coordinates this costs 0(q 2 ), rather than 
the O(qn) for (H|. 
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However, our experimental results show that DP-GLASSO is more efficient than GLASSO, so there 
are some other factors in play. When A is sparse, there are computational savings. If A has nq 
non-zeros, the cost per column reduces on average to O(Kq) from 0{q 2 ). For the formulation (|33|) A 
is ©11, which is sparse for large A. Hence for large A, GLASSO and dp-GLASSO have similar costs. 

For smaller values of A, the box-constrained QP (|37p is particularly attractive. Most of the 
coordinates in the optimal solution v will pile up at the boundary points {—A, A}, which means that 
the coordinates need not be updated frequently. For problem ([55)1 this number is also k, the number 
of non-zero coefficients in the corresponding column of the precision matrix. If k of the coordinates 
pile up at the boundary, then one full sweep of cyclical coordinate descent across all the coordinates 
will require updating gradients corresponding to the remaining q — k coordinates. Using similar 
calculations as before, this will cost 0(q(q — k)) operations per full cycle (since for small A, A will 
be dense). For the l\ regularized problem no such saving is achieved, and the cost is 0(q 2 ) per 
cycle. 

Note that to solve problem ([I]), we need to solve a QP of a particular type (|36|) or ([37)1 for a 
certain number of outer cycles (ie full sweeps across rows/columns). For every row/column update, 
the associated QP requires varying number of iterations to converge. It is hard to characterize all these 
factors and come up with precise estimates of convergence rates of the overall algorithm. However, 
we have observed that with warm-starts, on a relatively dense grid of As, the complexities given 
above are pretty much accurate for DP-GLASSO (with warmstarts) specially when one is interested 
in solutions with small / moderate accuracy. Our experimental results in Section 19.11 and Appendix 
Section [B] support our observation. 

We will now have a more critical look at the updates of the GLASSO algorithm and study their 
properties. 

7 GLASSO: Positive definiteness, Sparsity and Exact Inversion 

As noted earlier, GLASSO operates on W — it does not explicitly compute the inverse W _1 . It does 
however keep track of the estimates for 612 after every row/column update. The copy of retained 
by GLASSO along the row/column updates is not the exact inverse of the optimization variable W. 
Figure [2] illustrates this by plotting the squared-norm ||(0 — W _1 )||p as a function of the iteration 
index. Only upon (asymptotic) convergence, will be equal to W _1 . This can have important 
consequences. 



1 1 1 1 1 1 r~ 

100 200 300 
Iteration Index Iteration Index 



Figure 2: Figure illustrating some negative properties of GLASSO using a typical numerical example. [Left 
Panel] The precision matrix produced after every row/column update need not be the exact inverse of the 
working covariance matrix — the squared Frobenius norm of the error is being plotted across iterations. [Right 
Panel] The estimated precision matrix produced by GLASSO need not be positive definite along iterations; 
plot shows minimal eigen-value. 
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In many real- life problems one only needs an approximate solution to (fTj): 

• for computational reasons it might be impractical to obtain a solution of high accuracy; 

• from a statistical viewpoint it might be sufficient to obtain an approximate solution for that 
is both sparse and positive definite 

It turns out that the GLASSO algorithm is not suited to this purpose. 

Since the GLASSO is a block coordinate procedure on the covariance matrix, it maintains a positive 
definite covariance matrix at every row/column update. However, since the estimated precision matrix 
is not the exact inverse of W, it need not be positive definite. Although it is relatively straightforward 
to maintain an exact inverse of W along the row/column updates (via simple rank-one updates as 
before), this inverse W _1 need not be sparse. Arbitrary thresholding rules may be used to set 
some of the entries to zero, but that might destroy the positive-definiteness of the matrix. Since 
a principal motivation of solving ([1} is to obtain a sparse precision matrix (which is also positive 
definite), returning a dense W _1 to ([1} is not desirable. 

Figures [2] illustrates the above observations on a typical example. 

The dp-GLASSO algorithm operates on the primal (pj. Instead of optimizing the t\ regularized 
QP (|16j) . which requires computing 0^ , dp-GLASSO optimizes (|33|) . After every row/column update 
the precision matrix is positive definite. The working covariance matrix maintained by DP-GLASSO 
via W12 := S12 + 7 need not be the exact inverse of 0. Exact covariance matrix estimates, if required, 
can be obtained by tracking _1 via simple rank-one updates, as described earlier. 

Unlike GLASSO, DP-GLASSO (and P-GLASSO) return a sparse and positive definite precision matrix 
even if the row/column iterations arc terminated prematurely. 

8 Warm Starts and Path-seeking Strategies 

Since we seldom know in advance a good value of A, we often compute a sequence of solutions to ([T|) for 
a (typically) decreasing sequence of values Ai > A2 > • ■ • > \k- Warm-start or continuation methods 
use the solution at Ai as an initial guess for the solution at Aj+i, and often yield great efficiency. It 
turns out that for algorithms like GLASSO which operate on the dual problem, not all warm-starts 
necessarily lead to a convergent algorithm. We address this aspect in detail in this section. 

The following lemma states the conditions under which the row/column updates of the GLASSO 
algorithm will maintain positive definiteness of the covariance matrix W. 

Lemma 3. Suppose Z is used as a warm- start for the GLASSO algorithm. IfZ^O and ||Z— S||oo < A, 
then every row/column update of GLASSO maintains positive definiteness of the working covariance 
matrix W. 

Proof. Recall that the GLASSO solves the dual p9[) . Assume Z is partitioned as in and the pth 
row/column is being updated. Since Z >- 0, we have both 

Z11 >- and (z 22 - Z2i(Zii) _1 Zi 2 ) > 0. (38) 

Since Zn remains fixed, it suffices to show that after the row/column update, the expression (W22 — 
w 2 i(Zn) _1 Wi2) remains positive. Recall that, via standard optimality conditions we have W22 = 
S22 + A, which makes W22 > ^22 (since by assumption, I222 — S22 1 < A and Z22 > 0). Furthermore, 
W21 = S21 +7, where 7 is the optimal solution to the corresponding box-QP (|27|) . Since the starting 
solution Z21 satisfies the box-constraint l|27p i.e. ||z2i— S21 ||oo < A, the optimal solution of the QP (|2"7| 
improves the objective: 

W 2 l(Zii) _1 Wi2 < Z 2 l(Zii) _1 Zi2 

Combining the above along with the fact that W22 > ^22 we see 

W22 - w 2 i(Zn)~ 1 w 12 > 0, (39) 
which implies that the new covariance estimate W ^ 0. □ 
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Remark 3. If the condition ||Z — S||oo < A appearing in Lemma\3\ is violated, then the row/column 
update of GLASSO need not maintain PD of the covariance matrix W. 

We have encountered many counter-examples that show this to be true, see the discussion below. 

The R package implementation of GLASSO allows the user to specify a warm-start as a tuple 
(0o, Wo). This option is typically used in the construction of a path algorithm. 

If (©a, Wj) is provided as a warm-start for A' < A, then the GLASSO algorithm is not guaranteed to 
converge. It is easy to find numerical examples by choosing the gap A — A' to be large enough. Among 
the various examples we encountered, we briefly describe one here. Details of the experiment/data 
and other examples can be found in the online Appendix IA.1I We generated a data- matrix X„ xp , 
with n — 2,p = 5 with iid standard Gaussian entries. S is the sample covariance matrix. We 
solved problem ([I) using GLASSO for A = 0.9 x max^ 3 - |s,j|. We took the estimated covariance and 
precision matrices: Wj and ©a as a warm-start for the GLASSO algorithm with A' = A x 0.01. The 
GLASSO algorithm failed to converge with this warm-start. We note that |[Wa — S||oo = 0.0402 ^ A' 
(hence violating the sufficient condition in Lemma 0]) and after updating the first row/column via 
the GLASSO algorithm we observed that "covariance matrix" W has negative eigen-values — leading 
to a non-convergent algorithm. The above phenomenon is not surprising and easy to explain and 
generalize. Since Wa solves the dual (jTSJ), it is necessarily of the form Wj = S + T, for ||r||oo < A. 
In the light of Lemma [3J and also Remark [3J the warm-start needs to be dual-feasible in order to 
guarantee that the iterates W remain PD and hence for the sub-problems to be well defined convex 
programs. Clearly Wa does not satisfy the box-constraint ||Wa — < A', for A' < A. However, in 
practice the GLASSO algorithm is usually seen to converge (numerically) when A' is quite close to A. 

The following lemma establishes that any PD matrix can be taken as a warm-start for P-GLASSO 
or DP-GLASSOto ensure a convergent algorithm. 

Lemma 4. Suppose $ >- is a used as a warm-start for the P-GLASSO ( or DP-GLASSO ) algorithm. 
Then every row /column update of P-GLASSO (or DP-GLASSO ) maintains positive definiteness of the 
working precision matrix ©. 

Proof. Consider updating the pth row/column of the precision matrix. The condition $ >- is 
equivalent to both 

*n y and (022 - *2i(*n) _1 *i2) > 0. 

Note that the block <frn remains fixed; only the pth row/column of © changes. <p 2 \ gets updated to 
021; as docs &yi- From ([7} the updated diagonal entry # 2 2 satisfies: 

#22 - e2l(*ll) _1 012 = 7 ~~ TT > 0. 

(S 22 + A) 

Thus the updated matrix © remains PD. The result for the DP-GLASSO algorithm follows, since both 
the versions P-GLASSO and DP-GLASSO solve the same block coordinate problem. □ 

Remark 4. A simple consequence of Lemmas&and^is that the QPs arising in the process, namely 
the l\ regularized QPs U6\) and the box- constrained QPs and i3S\) are all valid convex 

programs, since all the respective matrices Wu, ©h an d ~^ix > ®u appearing in the quadratic 
forms are PD. 

As exhibited in Lemma HJ both the algorithms DP-GLASSO and P-GLASSO are guaranteed to 
converge from any positive-definite warm start. This is due to the unconstrained formulation of the 
primal problem (p}. 

GLASSO really only requires an initialization for W, since it constructs © on the fly. Likewise 
DP-GLASSO only requires an initialization for 0. Having the other half of the tuple assists in the 
block-updating algorithms. For example, GLASSO solves a series of lasso problems, where play 
the role as parameters. By supplying © along with W, the block-wise lasso problems can be given 
starting values close to the solutions. The same applies to DP-GLASSO. In neither case do the pairs 
have to be inverses of each other to serve this purpose. 
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If we wish to start with inverse pairs, and maintain such a relationship, we have described earlier 
how 0(p 2 ) updates after each block optimization can achieve this. One caveat for GLASSO is that 
starting with an inverse pair costs 0(p 3 ) operations, since we typically start with W = S + AI. For 
DP-GLASSO, we typically start with a diagonal matrix, which is trivial to invert. 



9 Experimental Results &; Timing Comparisons 

We compared the performances of algorithms GLASSO and dp-GLASSO (both with and without warm- 
starts) on different examples with varying (n,p) values. While most of the results are presented 
in this section, some are relegated to the online Appendix [BJ Section 19.11 describes some synthetic 
examples and Section [9.21 presents comparisons on a real-life micro-array data-set. 



9.1 Synthetic Experiments 

In this section we present examples generated from two different covariancc models — as characterized 
by the covariance matrix X or equivalently the precision matrix 0. We create a data matrix X nxp by 
drawing n independent samples from a p dimensional normal distribution MVN(0, S). The sample 
covariance matrix is taken as the input S to problem (TTJ). The two covariance models are described 
below: 

Type-1 The population concentration matrix = S _1 has uniform sparsity with approximately 77 
% of the entries zero. 

We created the covariance matrix as follows. We generated a matrix B with iid standard 
Gaussian entries, symmetrized it via i(B + B') and set approximately 77% of the entries of 
this matrix to zero, to obtain B (say). We added a scalar multiple of the p dimensional identity 
matrix to B to get the precision matrix = B + r)l pxp , with rj chosen such that the minimum 
eigen value of is one. 



Type-2 This example, taken from lYuan and Linl [2007j . is an auto-regressive process of order two 



the precision matrix being tri-diagonal: 



'0.5, if |j-i| = l, i = 2,...,(p-l); 

0. 25, if \j-i\ = 2, i = 3,...,(p- 2); 

1 , if i = j, i = 1, . . . , p; and 
otherwise. 



For each of the two set-ups Type-1 and Type-2 we consider twelve different combinations of (n,p): 

(a) p = 1000, n G {1500, 1000, 500}. 

(b) p = 800, n G {1000, 800, 500}. 

(c) p = 500, n £ {800, 500, 200}. 

(d) p = 200, n £ {500,200,50}. 

For every (n,p) we solved (JTJ on a grid of twenty A values linearly spaced in the log-scale, with 
Aj = 0.8* x {0.9A max }, i = 1, . . . , 20, where A max = max^ \sij\, is the off-diagonal entry of S with 
largest absolute value. A max is the smallest value of A for which the solution to ([T]) is a diagonal 
matrix. 

Since this article focuses on the GLASSO algorithm, its properties and alternatives that stem from 
the main idea of block-coordinate optimization, we present here the performances of the following 
algorithms: 



Dual-Cold GLASSO with initialization W = S + AI pxp , as suggested in iFriedman et al.l [2007 1. 
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Dual -War m The path-wise version of GLASSO with warm-starts, as suggested in iFriedman et al 



(2007 1 . Although this path- wise version need not converge in general, this was not a problem 



in our experiments, probably due to the fine-grid of A values. 
Primal-Cold dp-GLASSO with diagonal initialization = (diag(S) + AI) _1 . 
Primal- Warm The path-wise version of dp-GLASSO with warm-starts. 

We did not include P-GLASSO in the comparisons above since P-GLASSO requires additional matrix 
rank-one updates after every row/column update, which makes it more expensive. None of the above 
listed algorithms require matrix inversions (via rank one updates). Furthermore, dp-GLASSO and 
p-GLASSO arc quite similar as both are doing a block coordinate optimization on the dual. Hence 
we only included DP-GLASSO in our comparisons. We used our own implementation of the GLASSO 
and DP-GLASSO algorithm in R. The entire program is written in R, except the inner block-update 
solvers, which are the real work-horses: 



For G LASSO we used the lasso code crossProdLasso written in FORTRAN by IFriedman et al 
1200 



• For DP-GLASSO we wrote our own FORTRAN code to solve the box QP. 

An R package implementing DP-GLASSO will be made available in CRAN. 

In the figure and tables that follow below, for every algorithm, at a fixed A we report the total 
time taken by all the QPs — the l\ regularized QP for GLASSO and the box constrained QP for DP- 
GLASSO till convergence All computations were done on a Linux machine with model specs: Intcl(R) 
Xcon(R) CPU 5160 @ 3.00GHz. 

Convergence Criterion: Since DP-GLASSO operates on the the primal formulation and GLASSO 
operates on the dual — to make the convergence criteria comparable across examples we based it on 
the relative change in the primal objective values i.e. /(©) Q} across two successive iterations: 

/(Qfc ) o /(0 , fc - l) < TOL, (40) 

where one iteration refers to a full sweep across p rows /columns of the precision matrix (for DP-GLASSO 
) and covariance matrix (for GLASSO ); and TOL denotes the tolerance level or level of accuracy of 
the solution. To compute the primal objective value for the GLASSO algorithm, the precision matrix is 
computed from W via direct inversion (the time taken for inversion and objective value computation 
is not included in the timing comparisons). 

Computing the objective function is quite expensive relative to the computational cost of the 
iterations. In our experience convergence criteria based on a relative change in the precision matrix 
for DP-GLASSO and the covariance matrix for GLASSO seemed to be a practical choice for the examples 
we considered. However, for reasons we described above, we used criterion [40] in the experiments. 

Observations: Figure [4] presents the times taken by the algorithms to converge to an accuracy of 
TOL = 10 -4 on a grid of A values. 

The figure shows eight different scenarios with p > n, corresponding to the two different covariance 
models Type-1 (left panel) and Type-2 (right panel). It is quite evident that dp-GLASSO with warm- 
starts (Primal- Warm) outperforms all the other algorithms across all the different examples. All the 
algorithms converge quickly for large values of A (typically high sparsity) and become slower with 
decreasing A. For large p and small A, convergence is slow; however for p > n, the non-sparse end 
of the regularization path is really not that interesting from a statistical viewpoint. Warm-starts 
apparently do not always help in speeding up the convergence of GLASSO ; for example see Figure [4] 
with (n,p) = (500,1000) (Type 1) and (n,p) = (500,800) (Type 2). This probably further validates 
the fact that warm-starts in the case of GLASSO need to be carefully designed, in order for them 
to speed-up convergence. Note however, that GLASSO with the warm-starts prescribed is not even 
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guaranteed to converge — we however did not come across any such instance among the experiments 
presented in this section. 

Based on the suggestion of a referee we annotated the plots in Figure [4] with locations in the 
regularization path that are of interest. For each plot, two vertical dotted lines are drawn which 
correspond to the As at which the distance of the estimated precision matrix &\ from the population 
precision matrix is minimized wrt to the || • ||i norm (green) and || • \\p norm (blue). The optimal A 
corresponding to the || • ||i metric chooses sparser models than those chosen by \\-\\f', the performance 
gains achieved by DP-GLASSO seem to be more prominent for the latter A. 

Table Q] presents the timings for all the four algorithmic variants on the twelve different (n,p) 
combinations listed above for Type 1. For every example, we report the total time till convergence 
on a grid of twenty A values for two different tolerance levels: TOL 6 {10 -4 , 10 -5 }. Note that the 
DP-GLASSO returns positive definite and sparse precision matrices even if the algorithm is terminated 
at a relatively small/moderate accuracy level — this is not the case in GLASSO . The rightmost 
column presents the proportion of non-zeros averaged across the entire path of solutions ©a, where 
©a is obtained by solving {T]) to a high precision i.e. 10 -6 , by algorithms GLASSO and DP-GLASSO 
and averaging the results. 

Again we see that in all the examples dp-GLASSO with warm-starts is the clear winner among its 
competitors. For a fixed p, the total time to trace out the path generally decreases with increasing 
n. There is no clear winner between GLASSO with warm-starts and GLASSO without warm-starts. It 
is often seen that DP-GLASSO without warm-starts converges faster than both the variants of GLASSO 
(with and without warm-starts). 

Table [2] reports the timing comparisons for Type 2. Once again we see that in all the examples 
Primal- Warm turns out to be the clear winner. 

For n < p = 1000, we observe that Primal- Warm is generally faster for Type-2 than Type-1. 
This however, is reversed for smaller values of p £ {800,500}. Primal-Cold is has a smaller overall 
computation time for Type-1 over Type-2. In some cases (for example n < p = 1000), we see that 
Primal- Warm in Type-2 converges much faster than its competitors on a relative scale than in Type-1 

this difference is due to the variations in the structure of the covariance matrix. 



9.2 Micro-array Example 

We c o nsider the data-set int r oduce d in Alon et al. |l999| and further studied in IRothman et al 



2008j . iMazumder and Hastid |2012| . In this experiment, tissue samples were analyzed using an 



Affymetrix Oligonucleotide array. The data was processed, filtered and reduced to a subset of 2000 
gene expression values. The number of Colon Adenocarcinoma tissue samples is n = 62. For the 
purpose of the experiments presented in this section, we pre-screened the genes to a size of p = 725. 
We ob tained this subset of genes us ing the idea of exact covariance thresholding introduced in our 



paper [Mazumder and Hastid . |2012| . We thresholded the sample correlation matrix obtained from 



the 62 x 2000 microarray data-matrix into connected components with a threshold of 0.003643 
the genes belonging to the largest connected component formed our pre-screened gene pool of size 
p = 725. This (subset) data-matrix of size (n,p) = (62, 725) is used for our experiments. 

The results presented below in Table [3] show timing comparisons of the four different algorithms: 
Primal- Warm/Cold and Dual- Warm/Cold on a grid of fifteen A values in the log-scale. Once again 
we see that the Primal- Warm outperforms the others in terms of speed and accuracy. Dual- Warm 
performs quite well in this example. 



10 Conclusions 



This pap er explores some of th e apparent mysteries in the behavior of the GLASSO algorithm intro- 



duced in Friedman et al 



200?1 | . These have been explained by leveraging the fact that the GLASSO 
algorithm is solving the dual of the graphical lasso problem ([1]), by block coordinate ascent. Each 
block update, itself the solution to a convex program, is solved via its own dual, which is equivalent 



1 this is the largest value of the threshold for which the size of the largest connected component is smaller than 800 
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Type = 1 Type = 2 

n = 500/p= 1000 n = 500/p=1000 




0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 



n = 200/p = 500 n = 200/p = 500 




0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 




Figure 3: The timings in seconds for the four different algorithmic versions: GLASSO (with and without warm- 
starts) and dp-GLASSO (with and without warm-starts) for a grid of X values on the log-scale. [Left Panel] 
Covariance model for Type-1, [Right Panel] Covariance model for Type-2. The horizontal axis is indexed by 
the proportion of zeros in the solution. The vertical dashed lines correspond to the optimal X values for which 
the estimated errors \\&x — ©||i (green) and \\&x — &\\f (blue) are minimum. 
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p / n 


relative 
error (TOL) 


Total time (sees) to compute a path of solutions 
Dual- Cold Dual- Warm Primal- Cold Primal- Warm 


Average % 
Zeros in path 


1000 / 500 


10" 4 
IO" 5 


3550.71 
4706.22 


6592.63 
8835.59 


2558.83 
3234.97 


2005.25 
2832.15 


80.2 


1000 / 1000 


IO" 4 
10~ 5 


2788.30 
3597.21 


3158.71 
4232.92 


2206.95 
2710.34 


1347.05 
1865.57 


83.0 


1000 / 1500 


IO" 4 
IO" 5 


2447.19 
2764.23 


4505.02 
6426.49 


1813.61 
2199.53 


932.34 
1382.64 


85.6 


800 / 500 


10~ 4 
10~ 5 


1216.30 
1776.72 


2284.56 
3010.15 


928.37 
1173.76 


541.66 
798.93 


78.8 


800 / 800 


IO" 4 
IO" 5 


1135.73 
1481.36 


1049.16 
1397.25 


788.12 
986.19 


438.46 
614.98 


80.0 


800 / 1000 


IO" 4 
IO" 5 


1129.01 
1430.77 


1146.63 
1618.41 


786.02 
992.13 


453.06 
642.90 


80.2 


500 / 200 


io- 4 

IO" 5 


605.45 
811.58 


559.14 
795.43 


395.11 
520.98 


191.88 
282.65 


75.9 


500 / 500 


IO" 4 
IO" 5 


427.85 
551.11 


241.90 
315.86 


252.83 
319.89 


123.35 
182.81 


75.2 


500 / 800 


IO" 4 
10" 5 


359.78 
416.87 


279.67 
402.61 


207.28 
257.06 


111.92 
157.13 


80.9 


200 / 50 


10~ 4 
10" 5 


65.87 
92.04 


50.99 
75.06 


37.40 
45.88 


23.32 
35.81 


75.6 


200 / 200 


IO" 4 
IO" 5 


35.29 
45.90 


25.70 
33.23 


17.32 
22.41 


11.72 
17.16 


66.8 


200 / 300 


IO" 4 
10~ 5 


32.29 
38.37 


23.60 
33.95 


16.30 
20.12 


10.77 
15.12 


66.0 



Table 1: Table showing the performances of the four algorithms GLASSO (Dual-Warm/Cold) and dp-GLASSO 
(Primal-Warm /Cold) for the covariance model Type-1. We present the times (in seconds) required to compute 
a path of solutions to {JJ) (on a grid of twenty A values) for different (n,p) combinations and relative errors 
(as in 140V )- The rightmost column gives the averaged sparsity level across the grid of A values. dp-GLASSO 
with warm-starts is consistently the winner across all the examples. 
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p / n 


relative 
error (TOL) 


Total time (sees) to compute a path of solutions 
Dual- Cold Dual- Warm Primal- Cold Primal- Warm 


Average % 
Zeros in path 


1000 / 500 


10~ 4 


6093.11 


f A O O AO 

5483.03 


O A A f f* T 

3495.67 


1661.93 


75.6 


10~ 5 


7707.24 


7923.80 


A A A 1 O O 

4401.28 


2358.08 


1000 / 1000 


i A — 4 

10 


AT7Q HQ 

4/ <o.9o 


QCOO OQ 

dOOZ.ZS 


zo9 ( .oo 


11)15.84 


76.70 


in — 5 
10 


0054.21 


4714.80 


3444.79 


1593.54 


1000 / 1500 


1 O.— 4 
1U 


4 tav.Za 


ol ( 5.16 


nono on 
zo9o.o9 


1 a/it) nc 


78.5 


in — 5 
1U 


6171.01 


C A C O OA 

6958.29 


O /I o o o o 

3432.33 


1679.16 


800 / 500 


10~ 4 


2914.63 


3466.49 


1685.41 


1293.18 


74.3 


10~ 5 


O/ 1 T A TO 

3574.73 


/I C TO AT 

4572.97 


OAOO OA 

2083.20 


1893.22 


800 / 800 


in — 4 
10 


onoi KK 
ZUzl.Ou 


i nor; on 


11Q1 QC 

llol.oO 


DlO.UO 


74.4 


1U 


2521. Do 


opon /^o 

2o39.o2 


1 /I 1 C AC 

1415.95 


nor) no 

922.93 


800 / 1000 


1U 


oO ( 4.O0 


ZOOl.UD 


1oo4.oD 




75.9 


in — 5 
10 


4599.59 


ooro to 

3353.78 


O O G A r* o 

2260.58 


1353.28 


500 / 200 


10~ 4 


1200.24 


885.76 


718.75 


291.61 


70.5 


10~ 5 


1574.62 


1219.12 


876.45 


408.41 


500 / 500 


10~ 4 


575.53 


386.20 


323.30 


130.59 


72.2 


10~ 5 


730.54 


535.58 


421.91 


193.08 


500 / 800 


10" 4 


666.75 


474.12 


373.60 


115.75 


73.7 


10" 5 


852.54 


659.58 


485.47 


185.60 


200 / 50 


10~ 4 


110.18 


98.23 


48.98 


26.97 


73.0 


10" 5 


142.77 


133.67 


55.27 


33.95 


200 / 200 


10" 4 


50.63 


40.68 


23.94 


9.97 


63.7 


10" 5 


66.63 


56.71 


31.57 


14.70 


200 / 300 


10" 4 


47.63 


36.18 


21.24 


8.19 


65.0 


10~ 5 


60.98 


50.52 


27.41 


12.22 



Table 2: Table showing comparative timings of the four algorithmic variants of GLASSO and DP-GLASSO for 
the covariance model in Type-2. This table is similar to Table\T\ displaying results for Type-1. dp-GLASSO 
with warm-starts consistently outperforms all its competitors. 



relative Total time (sees) to compute a path of solutions 

error (TOL) Dual-Cold Dual- Warm Primal-Cold Primal- Warm 

HP 5 515.15 406.57 462.58 334.56 

10~ 4 976.16 677.76 709.83 521.44 



Table 3: Comparisons among algorithms for a microarray dataset with n — 62 and p = 725, for different 
tolerance levels (TOL). We took a grid of fifteen A values, the average % of zeros along the whole path is 90.8. 
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to a lasso problem. The optimization variable is W, the covariance matrix, rather than the target 
precision matrix 0. During the course of the iterations, a working version of is maintained, but it 
may not be positive definite, and its inverse is not W. Tight convergence is therefore essential, for 
the solution © to be a proper inverse covariance. There are issues using warm starts with GLASSO, 
when computing a path of solutions. Unless the sequence of As are sufficiently close, since the "warm 
start" s are not dual feasible, the algorithm can get into trouble. 

We have also developed two primal algorithms P-GLASSO and DP-GLASSO. The former is more 
expensive, since it maintains the relationship W = _1 at every step, an 0(p 3 ) operation per sweep 
across all row/columns. dp-GLASSO is similar in flavor to GLASSO except its optimization variable 
is 0. It also solves the dual problem when computing its block update, in this case a box-QP. This 
box-QP has attractive sparsity properties at both ends of the regularization path, as evidenced in 
some of our experiments. It maintains a positive definite © throughout its iterations, and can be 
started at any positive definite matrix. Our experiments show in addition that dp-GLASSO is faster 
than GLASSO. 

An R package implementing DP-GLASSO will be made available in CRAN. 
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A Online Appendix 

This section complements the examples provided in the paper with further experiments and illustra- 
tions. 

A.l Examples: Non-Convergence of GLASSO with warm-starts 

This section illustrates with examples that warm-starts for the GLASSO need not converge. This is a 
continuation of examples presented in Section [5] 

Example 1: 

We took (n,p) = (2, 5) and setting the seed of the random number generator in R as set.seed(2008) 
we generated a data-matrix X„ xp with iid standard Gaussian entries. The sample covariance matrix 
S is given below: 



/ 0.03597652 
0.03597652 
0.10585853 
-0.08360659 
\ 0.13667246 



0.03792221 
0.03792221 
0.11158361 
-0.08812823 
0.14406402 



0.1058585 
0.1058585 
0.3114818 
-0.2460069 
0.4021497 



-0.08360659 
-0.08360659 
-0.24600689 
0.19429514 
-0.31761603 



0.1366725 \ 
0.1366725 
0.4021497 
-0.3176160 
0.5192098 / 



With q denoting the maximum off-diagonal entry of S (in absolute value), we solved ([T]) using 
GLASSO at A = 0.9 x q. The covariance matrix for this A was taken as a warm-start for the GLASSO 
algorithm with A' = A x 0.01. The smallest eigen- value of the working covariance matrix W produced 
by the GLASSO algorithm, upon updating the first row/column was: —0.002896128, which is clearly 
undesirable for the convergence of the algorithm GLASSO . This is why the algorithm GLASSO breaks 
down. 



Example 2: 

The example is similar to above, with (n,p) 



(10,50), the seed of random number generator in R 



being set to set.seed(2008) and X„ X p is the data-matrix with iid Gaussian entries. If the covariance 
matrix Wa which solves problem |T]) with A = 0.9 x maxj^j \sij\ is taken as a warm-start to the 
GLASSO algorithm with A' = A x 0.1 — the algorithm fails to converge. Like the previous example, 
after the first row/column update, the working covariance matrix has negative eigen-values. 



B Further Experiments and Numerical Studies 

This section is a continuation to Section [HI in that it provides further examples comparing the 
performance of algorithms GLASSO and dp-GLASSO . The experimental data is generated as follows. 
For a fixed value of p, we generate a matrix A pxp with random Gaussian entries. The matrix is 
symmetrized by A <— (A + A')/2. Approximately half of the off-diagonal entries of the matrix are 
set to zero, uniformly at random. All the eigen-values of the matrix A are lifted so that the smallest 
eigen- value is zero. The noiseless version of the precision matrix is given by © = A + Tl pX p- We 
generated the sample covariance matrix S by adding symmetric positive semi-definite random noise 
N to i.e. S = _1 + N, where this noise is generated in the same manner as A. We considered 

four different values of p e {300, 500, 800, 1000} and two different values of r G {1,4}. 

For every p, r combination we considered a path of twenty A values on the geometric scale. For 
every such case four experiments were performed: Primal-Cold, Primal- Warm, Dual-Cold and Dual- 
Warm (as described in Section [5]). Each combination was run 5 times, and the results averaged, to 
avoid dependencies on machine loads. Figure [4] shows the results. Overall, dp-GLASSO with warm 
starts performs the best, especially at the extremes of the path. We gave some explanation for this 
in Section [6l For the largest problems (p = 1000) their performances are comparable in the central 
part of the path (though DP-GLASSO dominates), but at the extremes dp-GLASSO dominates by a 
large margin. 
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Figure 4: The timings in seconds for the four different algorithmic versions GLASSO (with and without 
warm-starts) and DP-GLASSO (with and without warm-starts) for a grid of twenty A values on the 
log-scale. The horizontal axis is indexed by the proportion of zeros in the solution. 



21 



